library(foreign)
library(ggplot2)
library(ggeffects)
library(ggpubr)
library(gridExtra)
library(arm)
library(Rmisc)
library(RColorBrewer)
library(car)
library(lmtest)
library(sandwich)
library(cregg)
library(dplyr)



# Conjoint experiments
load("Conjoint_data.rda")

## Figure 1
amces <- cj(Conjoint_data, chosen2 ~ `National economy` + `Local economy` + `Fiscal burden` + `Corporate culture` + `National culture`, id=~IDs, alpha=0.05)
plot(amces, xlab="Change in Pr(Preferred future)") + theme(legend.position="none") + theme(axis.text=element_text(size=10))  + scale_colour_manual(values=rep("black", 5)) 


## Figure A4
amces_f <- cj(Conjoint_data, chosen2 ~ `Security` + `National economy` + `Local economy` + `Fiscal burden` + `Corporate culture` + `National culture`, id=~IDs, alpha=0.05)
plot(amces_f, xlab="Change in Pr(Preferred future)") + theme(legend.position="none") + theme(axis.text=element_text(size=10))  + scale_colour_manual(values=rep("black", 6)) 


## Figure A5
Conjoint_data$college_above_f[Conjoint_data$college_above == 1] <- "College or above"
Conjoint_data$college_above_f[Conjoint_data$college_above == 0] <- "Below college"
Conjoint_data$college_above_f <- factor(Conjoint_data$college_above_f)

mm_by_edu <- cj(Conjoint_data, chosen2 ~ `National economy` + `Local economy` + `Fiscal burden` + `Corporate culture` + `National culture`, id=~IDs, estimate="mm_differences",by=~college_above_f)
plot(mm_by_edu, group="college_above_f", vline=0, legend_title=NULL) + theme(legend.position="none") + theme(legend.position="none") + theme(axis.text=element_text(size=10)) + scale_colour_manual(values=rep("black", 5))


## Figure A6
Conjoint_data$age[Conjoint_data$age_normal >= 0.49] <- "Old"
Conjoint_data$age[Conjoint_data$age_normal < 0.49] <- "Young"
Conjoint_data$age <- factor(Conjoint_data$age)

mm_by_old <- cj(Conjoint_data, chosen2 ~ `National economy` + `Local economy` + `Fiscal burden` + `Corporate culture` + `National culture`, id=~IDs, estimate="mm_differences",by=~age)
plot(mm_by_old, group="age", vline=0, legend_title=NULL) + theme(legend.position="none") + theme(legend.position="none") + theme(axis.text=element_text(size=10))  + scale_colour_manual(values=rep("black", 5)) 


## Figure A7
Conjoint_data$income[Conjoint_data$income_normal >= 0.2857] <- "Rich"
Conjoint_data$income[Conjoint_data$income_normal < 0.2857] <- "Poor"
Conjoint_data$income <- factor(Conjoint_data$income)

mm_by_rich <- cj(Conjoint_data, chosen2 ~ `National economy` + `Local economy` + `Fiscal burden` + `Corporate culture` + `National culture`, id=~IDs, estimate="mm_differences",by=~income)
plot(mm_by_rich, group="income", vline=0, legend_title=NULL) + theme(legend.position="none") + theme(legend.position="none") + theme(axis.text=element_text(size=10)) + scale_colour_manual(values=rep("black", 5)) 




# Vignette experiments (2016)
load("Vignette_data.rda")


## Figure 2
mean <- summarySE(Vignette_data, measurevar="opposition_w2_r", groupvars=c("Treatment_w2_s"), na.rm=TRUE, conf.interval=0.95)
mean

p = ggplot(mean, aes(y=opposition_w2_r, x=Treatment_w2_s, colour=Treatment_w2_s)) + geom_hline(yintercept=3, linetype="dashed") + geom_errorbar(aes(ymin=opposition_w2_r-ci, ymax=opposition_w2_r+ci), width=.1, position=position_dodge(.5), color="black") + geom_point(aes(shape=Treatment_w2_s), size=3,color="black") + scale_shape_manual(values=c(0,1,16,5,2,17),name="Treatment groups", labels=c("Control","Economic (positive)","Economic (negative)","Fiscal (negative)","Culture (positive)","Culture (negative)")) + scale_y_continuous(limits=c(1, 4)) + scale_x_discrete(labels=c("","","","","","","","","","","")) + labs(title="Mean support rates") + ylab("Mean scores") + xlab("") + theme_classic() + theme(legend.key=element_rect(colour=NA))


## Figure 3
Vignette_data$college_above_f[Vignette_data$college_above == 1] <- "College or above"
Vignette_data$college_above_f[Vignette_data$college_above == 0] <- "Below college"
Vignette_data$college_above_f <- factor(Vignette_data $college_above_f)

reg1 <- lm(opposition_w2_r ~ college_above_f*Treatment_w2_s + age_normal + male + income_normal + unemployment_real + LDP  + for_ami_friend_d + RWA_factor + ethno_nor, data=Vignette_data)
ggpredict(reg1,c("Treatment_w2_s","college_above_f")) %>% plot(colors="bw", dot.size=3.5) + geom_hline(yintercept=3, linetype="dashed") + labs(x = "", y = "Treatment effects by respondents' educational attainment", title = "") + theme_classic() + theme(legend.key=element_rect(colour=NA)) + theme(legend.title=element_blank(),legend.text=element_text(size=12)) + theme(axis.text.x=element_text(angle=45,hjust=1,size=13)) + scale_y_continuous(limits=c(1,4))
summary(reg1)


## Figure 4
quantile(Vignette_data$age_normal, na.rm=TRUE)
Vignette_data$age[Vignette_data$age_normal >= 0.52] <- "Old"
Vignette_data$age[Vignette_data$age_normal < 0.52] <- "Young"
Vignette_data$age <- as.factor(Vignette_data$age)

reg2 <- lm(opposition_w2_r ~ age*Treatment_w2_s + male + income_normal + college_above + unemployment_real + LDP  + for_ami_friend_d + RWA_factor + ethno_nor, data=Vignette_data)
ggpredict(reg2,c("Treatment_w2_s","age")) %>% plot(colors="bw", dot.size=3.5) + geom_hline(yintercept=3, linetype="dashed") + labs(x = "", y = "Treatment effects by respondents' age", title = "") + theme_classic() + theme(legend.key=element_rect(colour=NA)) + theme(legend.title=element_blank(),legend.text=element_text(size=12)) + theme(axis.text.x=element_text(angle=45,hjust=1,size=13)) + scale_y_continuous(limits=c(1,4)) 
summary(reg2)


## Figure 5
quantile(Vignette_data$income_normal, na.rm=TRUE)
Vignette_data$income[Vignette_data$income_normal >= 0.2857] <- "Rich"
Vignette_data$income[Vignette_data$income_normal < 0.2857] <- "Poor"
Vignette_data$income <- as.factor(Vignette_data$income)

reg3 <- lm(opposition_w2_r ~ income*Treatment_w2_s + male + age_normal + college_above + unemployment_real + LDP  + for_ami_friend_d + RWA_factor + ethno_nor, data=Vignette_data)
ggpredict(reg3,c("Treatment_w2_s","income")) %>% plot(colors="bw", dot.size=3.5) + geom_hline(yintercept=3, linetype="dashed") + labs(x = NULL, y = "Treatment effects by respondents' income", title = "") + theme_classic() + theme(legend.key=element_rect(colour=NA)) + theme(legend.title=element_blank(),legend.text=element_text(size=12)) + theme(axis.text.x=element_text(angle=45,hjust=1,size=13)) + scale_y_continuous(limits=c(1,4))
summary(reg3)




# Vignette experiments (2020)
load("Vignette_data2020.rda")


## Figure 6
mean2020 <- summarySE(Vignette_data2020, measurevar="support_r", groupvars=c("treatment"), na.rm=TRUE, conf.interval=0.95)
mean2020

p_2020 = ggplot(mean2020, aes(y=support_r, x=treatment, colour=treatment)) + geom_hline(yintercept=3, linetype="dashed") + geom_errorbar(aes(ymin=support_r-ci, ymax=support_r +ci), width=.1, position=position_dodge(.5), color="black") + geom_point(aes(shape=treatment), size=3,color="black") + scale_shape_manual(values=c(0,1,16),name="Treatment groups", labels=c("Control","Fiscal (negative)","Fiscal (positive)")) + scale_y_continuous(limits=c(1, 4)) + scale_x_discrete(labels=c("","","")) + labs(title="Mean support rates") + ylab("Mean scores") + xlab("") + theme_classic() + theme(legend.key=element_rect(colour=NA)) 

Vignette_data2020$employed[Vignette_data2020$employment == 1]  <- "Employed"
Vignette_data2020$employed[Vignette_data2020$employment >= 2 | Vignette_data2020$employment >= 7 ]  <- "Unemployed"
Vignette_data2020$employed <- factor(Vignette_data2020$employed)

Vignette_data2020$college_above_f[Vignette_data2020$college == 1] <- "College or above"
Vignette_data2020$college_above_f[Vignette_data2020$college == 0] <- "Below college"
Vignette_data2020$college_above_f <- factor(Vignette_data2020$college_above_f)

reg2020_college <- lm(support_r ~ college_above_f*treatment + age + women + income + employed, data=Vignette_data2020)
ggpredict(reg2020_college,c("treatment","college_above_f")) %>% plot(colors="bw", dot.size=3.5) + geom_hline(yintercept=3, linetype="dashed") + labs(x = "", y = "Treatment effects by respondents' educational attainment", title = "") + theme_classic() + theme(legend.key=element_rect(colour=NA)) + theme(legend.title=element_blank(),legend.text=element_text(size=12)) + theme(axis.text.x=element_text(angle=45,hjust=1,size=13)) + scale_y_continuous(limits=c(1,4))


median(Vignette_data2020$age, na.rm=TRUE)
Vignette_data2020$old_young[Vignette_data2020$age >= 52] <- "Old"
Vignette_data2020$old_young[Vignette_data2020$age < 52] <- "Young"
Vignette_data2020$old_young <- factor(Vignette_data2020$old_young)

reg2020_age <- lm(support_r ~ old_young*treatment + college_above_f + women + income + employed, data=Vignette_data2020)
ggpredict(reg2020_age,c("treatment","old_young")) %>% plot(colors="bw", dot.size=3.5) + geom_hline(yintercept=3, linetype="dashed") + labs(x = "", y = "Treatment effects by respondents' age", title = "") + theme_classic() + theme(legend.key=element_rect(colour=NA)) + theme(legend.title=element_blank(),legend.text=element_text(size=12)) + theme(axis.text.x=element_text(angle=45,hjust=1,size=13)) + scale_y_continuous(limits=c(1,4))


median(Vignette_data2020$income, na.rm=TRUE)
Vignette_data2020$rich_poor[Vignette_data2020$income >= 3] <- "Rich"
Vignette_data2020$rich_poor[Vignette_data2020$income < 3] <- "Poor"
Vignette_data2020$rich_poor <- factor(Vignette_data2020$rich_poor)

reg2020_income <- lm(support_r ~ rich_poor*treatment + age + college_above_f + women + employed, data=Vignette_data2020)
ggpredict(reg2020_income,c("treatment","rich_poor")) %>% plot(colors="bw", dot.size=3.5) + geom_hline(yintercept=3, linetype="dashed") + labs(x = "", y = "Treatment effects by respondents' income", title = "") + theme_classic() + theme(legend.key=element_rect(colour=NA)) + theme(legend.title=element_blank(),legend.text=element_text(size=12)) + theme(axis.text.x=element_text(angle=45,hjust=1,size=13)) + scale_y_continuous(limits=c(1,4))



